home *** CD-ROM | disk | FTP | other *** search
/ Languguage OS 2 / Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO / language / dino / dino_bot.1 / examples / complex / block.d next >
Encoding:
Text File  |  1991-04-15  |  13.5 KB  |  510 lines
  1. /* Copyright, 1990, Regents of the University of Colorado */
  2. /* This program is a fairly long example, designed to show that complex programs
  3.  * can be written in DINO. */
  4.  
  5. /****************************************************************************
  6.  * Solves system of linear equations with a block bordered structure
  7.  * (See Rosin, Schnabel, and Weaver's "Expressing Complex Parallel
  8.  * Algorithms in Dino" p. 2).
  9.  *
  10.  * General Form:
  11.  *
  12.  *  A1          B       x1       f1
  13.  *     A2       B       x2       f2
  14.  *        ..    ..  *   ..   =   ..
  15.  *           Aq Bq      xq       fq
  16.  *  C  C  .. Cq P       xqp      fqp
  17.  *
  18.  * Where:
  19.  *
  20.  *  A   is nxn
  21.  *  B   is nxm
  22.  *  C   is mxn
  23.  *  P   is mxm
  24.  *  x   is nx1
  25.  *  xqp is mx1
  26.  *  f   is nx1
  27.  *  fqp is mx1
  28.  *
  29.  * Command Line:
  30.  *
  31.  *  block <filename>
  32.  *
  33.  * where <filename> is the name of the file describing the linear system
  34.  * as specified in function read_data.
  35.  ***************************************************************************/
  36.  
  37. #include <stdio.h>
  38. #include <math.h>
  39. #include <string.h>             /* ==> Standard include files can be used
  40.                                        in the normal manner. */
  41.  
  42. #include "dino.h"               /* ==> The directory containing dino.h is
  43.                                        automatically searched. */
  44.  
  45. map byRow     = [block][compress];
  46. map byBlock   = [block];
  47. map byElement = [block];
  48. map byCol     = [compress][block];
  49.  
  50. #define N 2
  51. #define M 3
  52. #define Q 2
  53.  
  54. environment solve[M:id] {
  55.  
  56.     composite form_J( in sumCW, P )
  57.     double distributed sumCW[M][M] map byRow;
  58.     double distributed P[M][M] map byRow;
  59.     {
  60.         int i;
  61.         for( i=0; i<M; i++ )
  62.             P[id][i] = P[id][i] - sumCW[id][i];
  63.     };
  64.  
  65.     composite form_b( b, in sumCz, in fqp )
  66.     double distributed b[M]     map byElement;
  67.     double distributed sumCz[M] map byElement;
  68.     double distributed fqp[M]   map byElement;
  69.     {
  70.         b[id] = fqp[id] - sumCz[id];
  71.     };
  72.  
  73.     composite dist_lu( A, RowPerm )
  74.     double distributed A[M][M] map byCol;
  75.     int distributed RowPerm[M] map all;
  76.     {
  77.         double distributed Mult[M] map all;
  78.         int i, j, k, pivrow;
  79.         double piv, temp;
  80.  
  81.         RowPerm[M-1] = M-1;
  82.         for( k=0; k<M-1; k++ ) {
  83.             if( k == id ) {
  84.                 /* select pivot and swap in pivot column */
  85.                 piv = A[k][k];
  86.                 pivrow = k;
  87.  
  88.                 /* find largest row element in kth column (partial pivot) */
  89.                 for( i=k+1; i<M; i++ )
  90.                     if( fabs(A[i][k]) > piv ) {
  91.                         piv = A[i][k];
  92.                         pivrow = i;
  93.                     }
  94.  
  95.                 /* swap row position in kth column */
  96.                 if( pivrow != k ) {
  97.                     temp = A[k][k];
  98.                     A[k][k] = A[pivrow][k];
  99.                     A[pivrow][k] = temp;
  100.                 }
  101.  
  102.                 /* calculate l (multipliers) and row swap info. */
  103.                 RowPerm[k] = pivrow;
  104.                 for (i=k+1; i<M; i++)
  105.                     Mult[i] = A[i][k] /= -A[k][k];
  106.  
  107.                 /* broadcast Mult and RowPerm */
  108.                 Mult[<k+1, M-1>]#{ solve[] - solve[id] } =  Mult[<k+1, M-1>];
  109.                 RowPerm[k]#{ solve[] - solve[id] } = RowPerm[k];
  110.             }
  111.  
  112.             else {
  113.                 /* receive broadcast here */
  114.                 Mult[<k+1, M-1>] = Mult[<k+1, M-1>]#{ solve[k] };
  115.                 RowPerm[k] = RowPerm[k]#{ solve[k] };
  116.             }
  117.  
  118.             /* eliminate in your column greater than k */
  119.             pivrow = RowPerm[k];
  120.             if( id>k ) {
  121.                 if( pivrow != k ) {
  122.                     temp = A[k][id];
  123.                     A[k][id] = A[pivrow][id];
  124.                     A[pivrow][id] = temp;
  125.                 }
  126.                 for( i=k+1; i<M; i++ )
  127.                     A[i][id] += Mult[i]*A[k][id]; /* was k */
  128.             }
  129.             /* swap factors in lower triangular */
  130.             else if( id<k )
  131.                 for( i=k; i<M; i++ ) {
  132.                     temp = A[i][id];
  133.                     A[i][id] = A[pivrow][id];
  134.                     A[pivrow][id] = temp;
  135.                 }
  136.         }
  137.     };
  138.  
  139.     composite dist_solve( in lu, out x, in b, in RowPerm )
  140.     double distributed lu[M][M]     map byRow;
  141.     double distributed x[M]         map byElement;
  142.     double distributed b[M]         map all;
  143.     int distributed RowPerm[M]      map all;
  144.     {
  145.         double distributed y[M]     map byElement;
  146.         double temp;
  147.         int i, swap;
  148.  
  149.         /* permute b according to RowPerm */
  150.         for( i=0; i<=M-2; i++ ) {
  151.             swap = RowPerm[i];
  152.             if( swap != i ) {
  153.                 temp = b[i];
  154.                 b[i] = b[swap];
  155.                 b[swap] = temp;
  156.             }
  157.         }
  158.  
  159.         /* solve ly = b */
  160.         y[id] = b[id];
  161.         for( i=0; i<=id-1; i++ )
  162.             y[id] += lu[id][i]*y[i]#{solve[i]};
  163.         y[id]#{solve[<id+1,M-1>]} = y[id];
  164.  
  165.         /* solve ux=y */
  166.         x[id] = y[id];
  167.         for( i=M-1; i>id; i-- )
  168.             x[id] -= lu[id][i]*x[i]#{solve[i]};
  169.         x[id] = x[id] / lu[id][id];
  170.         x[id]#{solve[<0,id-1>]} = x[id];
  171.     };
  172. }
  173.  
  174. environment node[Q:id] {
  175.  
  176.     map slice = [block][compress][compress];
  177.  
  178.     double distributed W[Q][N][M] map slice;  /* = A^-1 B*/
  179.     double distributed z[Q*N]     map byBlock;
  180.                                 /* ==> Distributed variables may be declared
  181.                                        inside an environment, in which case
  182.                                        they are accessable by all functions and
  183.                                        composite procedures within that
  184.                                        environment. */
  185.  
  186.     /*
  187.      * negate v of length l
  188.      */
  189.     neg_vec( v,l )
  190.     double v[];
  191.     int l;
  192.     {
  193.         int i;
  194.         for( i=0; i<=l-1; i++ )
  195.             v[i] = v[i] * -1.0;
  196.     };
  197.  
  198.     /*
  199.      * v will contain (v - x) where v & x are 1 x l
  200.      */
  201.     sub_vec( v,x,l ) /* v <= v - x, l is length of vectors */
  202.     double v[];
  203.     double x[];
  204.     int l;
  205.     {
  206.         int i;
  207.         for( i=0; i<=l-1; i++ )
  208.             v[i] = v[i] - x[i];
  209.     };
  210.  
  211.     /*
  212.      * v will contain (v - x) where v & x are 1 x l
  213.      */
  214.     add_vec( v,x,l ) /* v <= v - x, l is length of vectors */
  215.     double v[];
  216.     double x[];
  217.     int l;
  218.     {
  219.         int i;
  220.         for( i=0; i<=l-1; i++ )
  221.             v[i] = v[i] + x[i];
  222.     };
  223.  
  224.     /*
  225.      * returns (A*B) where A is mxn and B nxm and T is mxm
  226.      */
  227.     mat_mat_mult(A,B,T)
  228.     double A[M][N], B[N][M], T[M][M];
  229.     {
  230.         int row, col, i;
  231.         for( row=0; row<=M-1; row++ )
  232.             for( col=0; col<=M-1; col++ ) {
  233.                 T[row][col] = 0;
  234.                 /* inner product */
  235.                 for( i=0; i<=N-1; i++ )
  236.                     T[row][col] += A[row][i]*B[i][col];
  237.             }
  238.     };
  239.  
  240.     /*
  241.      * Does lu decomposition on matrix A
  242.      * Note: o uses partial pivoting
  243.      *       o assume A is nxn since seq_lu is only used on A to find z
  244.      *         (see above mentioned paper)
  245.      *       o after call to seq_lu A will contain lu decomposition
  246.      */
  247.     void seq_lu( A, p )
  248.     double A[N][N];
  249.     int *p;
  250.     {
  251.         double temp, mult;
  252.         int i, diag, col, pivot;
  253.         for( diag=0; diag<=N-2; diag++ ) {
  254.  
  255.             /* look for largest leading row element, i.e. pivot */
  256.             pivot = diag;
  257.             for( i=diag; i<=N-1; i++ ) {
  258.                 if( A[i][diag] > A[pivot][diag] )
  259.                     pivot = i;
  260.             }
  261.             p[diag] = pivot;
  262.  
  263.             /* new pivot found? */
  264.             if( pivot != diag ) {
  265.                 /* information about row swaps stored in p such that
  266.                    p[diag] = pivot means that row pivot and diag were
  267.                    swapped at diagth iteration */
  268.                 for( i=0; i<=N-1; i++ ) { /* swap rows */
  269.                     temp = A[diag][i];
  270.                     A[diag][i] = A[pivot][i];
  271.                     A[pivot][i] = temp;
  272.                 }
  273.             }
  274.  
  275.             /* elimate leading columns */
  276.             for( i=diag+1; i<=N-1; i++ ) {
  277.                 /* compute l */
  278.                 mult = A[i][diag] /= -A[diag][diag];
  279.                 for( col=diag+1; col<=N-1; col++ )
  280.                     A[i][col] += mult*A[diag][col];
  281.             }
  282.         }
  283.         p[N-1] = N-1; /* last row can't be swapped at last iteration */
  284.     };
  285.  
  286.     /*
  287.      * solves: lu x = b
  288.      * Note: o lu is nxn (see seq_lu) and naturally x, b are 1xn
  289.      *   o result found in x after call to seq_solve
  290.      */
  291.     void seq_solve( lu, x, b, p ) /* solves lu x = b */
  292.     double lu[N][N];
  293.     double x[N];
  294.     double b[N];
  295.     int p[N];
  296.     {
  297.         int i, col, diag;
  298.  
  299.         double y[N];
  300.         double temp;
  301.  
  302.         /* solve ly = b */
  303.         for( diag=0; diag<=N-2; diag++ ) {
  304.             if( p[diag] != diag ) {
  305.                 temp = b[diag];
  306.                 b[diag] = b[p[diag]];
  307.                 b[p[diag]] = temp;
  308.             }
  309.         }
  310.  
  311.         y[0] = b[0];
  312.         for( diag=1; diag<=N-1; diag++ ) {
  313.             y[diag] = b[diag];
  314.             for( col=0; col<=diag-1; col++ )
  315.                 y[diag] += lu[diag][col]*y[col];
  316.         }
  317.  
  318.         /* solve ux=y */
  319.         x[N-1] = y[N-1]/lu[N-1][N-1];
  320.         for( diag=N-2; diag>=0; diag-- ) {
  321.             x[diag] = y[diag];
  322.             for( col=diag+1; col<=N-1; col++ )
  323.                 x[diag] -= lu[diag][col]*x[col];
  324.             x[diag] /= lu[diag][diag];
  325.         }
  326.  
  327.     };
  328.  
  329.     composite factor_A(in A, in f, in B)
  330.     double distributed A[Q][N][N] map slice;
  331.     double distributed f[Q*N]     map byBlock;
  332.     double distributed B[Q][N][M]  map slice;
  333.     {
  334.         int i,r;
  335.         int p[N];
  336.         double tempW[N];
  337.         double tempB[N];
  338.  
  339.         /* decompose each block */
  340.         seq_lu(A[id], p);
  341.         seq_solve(A[id], &z[id*N], &f[id*N], p);   /* z = A inv f */
  342.  
  343.         for (i=0; i<M; i++) {                      /* W = A inv B */
  344.             tempB[] = B[id][][i];
  345.             seq_solve(A[id], tempW, tempB, p);
  346.             W[id][][i] = tempW[];
  347.         }
  348.     };
  349.  
  350.     composite form_sums( in C, out sumCW, out sumCz )
  351.     double distributed C[Q][M][N]     map slice;
  352.     double distributed sumCW[M][M]    map byRow;
  353.     double distributed sumCz[M]       map byElement;
  354.     {
  355.         double T[M];
  356.         double T2[M];
  357.         double Temp[M][M];
  358.         int i;
  359.         int r,c;
  360.  
  361.         /* sum C*W  = J*/
  362.         mat_mat_mult(C[id], W[id], Temp);
  363.  
  364.         Temp[][] = gsum( Temp[][] )#;
  365.  
  366.         for( i=id*M/Q; i<=id*M/Q + M/Q; i++) {
  367.             sumCW[i][] = Temp[i][];
  368.         }
  369.  
  370.         if( id==Q-1 ) {
  371.             sumCW[M%Q+1][] = Temp[M%Q+1][];
  372.         }
  373.  
  374.         /* sum c*z */
  375.         for( r=0; r<=M-1; r++ ) {
  376.             T[r] = 0;
  377.             for( c=0; c<=N-1; c++ )
  378.                 T[r] += C[id][r][c]* z[id*N+c];
  379.         }
  380.  
  381.         T2[] = gsum(T[])#;
  382.  
  383.         for( i=id*M/Q; i<=id*M/Q + M/Q; i++)
  384.             sumCz[i] = T2[i];
  385.  
  386.         if( id==Q-1 )
  387.             sumCz[M%Q+1] = T2[M%Q+1];
  388.     };
  389.  
  390.     composite compute_xi( in xqp, out x )
  391.     double distributed xqp[M]     map all;
  392.     double distributed x[Q*N]     map byBlock;
  393.     {
  394.         int r,c;
  395.  
  396.         for( r=0; r<=N-1; r++ ) {
  397.             x[id*N+r] = 0;
  398.             for( c=0; c<=M-1; c++ )
  399.                 x[id*N+r] += W[id][r][c]*xqp[c];
  400.         }
  401.         neg_vec( &x[id*N], N );
  402.         add_vec( &x[id*N], &z[id*N], N );
  403.     };
  404. }
  405.  
  406. environment host{
  407.     double A[Q][N][N], B[Q][N][M], C[Q][M][N], P[M][M],
  408.            fqp[M], f[Q*N], xqp[M], x[Q*N], b[M], RowPerm[M], sumCW[M][M],
  409.            sumCz[M];
  410.     FILE *InFile;
  411.     char filename[256];
  412.  
  413.     void read_dvector( InFile, v, size )
  414.     FILE *InFile;
  415.     double *v;
  416.     int size;
  417.     {
  418.         int i;
  419.         double temp;
  420.  
  421.         for( i=0; i<=size-1; i++ )
  422.             fscanf( InFile, "%lf\n", &v[i] );
  423.     };
  424.  
  425.     void read_dmatrix( InFile, m, row, col )
  426.     FILE *InFile;
  427.     double *m;
  428.     int row, col;
  429.     {
  430.         int r, c;
  431.         for( r=0; r<=row-1; r++ )
  432.             read_dvector( InFile, (double *) (m+r*col), col );
  433.     };
  434.  
  435.     void read_dmatrices( InFile, ms, size, row, col )
  436.     FILE *InFile;
  437.     double *ms;
  438.     int size, row, col;
  439.     {
  440.         int i;
  441.  
  442.         for( i=0; i<=size-1; i++ )
  443.             read_dmatrix( InFile, (double *) (ms+i*row*col), row, col );
  444.     };
  445.  
  446.     /*
  447.      * read_data reads an ascii file containing the linear system.
  448.      * File must be in the following form (n, m, q must be set in program):
  449.      *
  450.      *  A1 <CR>
  451.      *  ...
  452.      *  Aq <CR>
  453.      *  B1 <CR>
  454.      *  ...
  455.      *  Bq <CR>
  456.      *  C1 <CR>
  457.      *  ...
  458.      *  Cq <CR>
  459.      *  P <CR>
  460.      *  f <CR>
  461.      *  fqp <CR>
  462.      *
  463.      * Note:
  464.      *
  465.      *  Matrices must be listed in row order.
  466.      *
  467.      */
  468.     void read_data( )
  469.     void;
  470.     {
  471.         InFile = fopen( filename, "r" );
  472.         read_dmatrices( InFile, A, Q, N, N );
  473.         read_dmatrices( InFile, B, Q, N, M );
  474.         read_dmatrices( InFile, C, Q, M, N );
  475.         read_dmatrix( InFile, P, M, M );
  476.         read_dvector( InFile, f, Q*N );
  477.         read_dvector( InFile, fqp, M );
  478.     };
  479.  
  480.     print_results()
  481.     {
  482.         int i, ii;
  483.  
  484.         printf( "Solution to System: \n" );
  485.         /* print x1 to xq */
  486.         for( i=0; i<=(Q*N)-1; i++ )
  487.             printf( "%lf\n", x[i] );
  488.         /* print xqp */
  489.         for( i=0; i<=M-1; i++ )
  490.             printf( "%lf\n", xqp[i] );
  491.     };
  492.  
  493.     void main( argc, argv )
  494.     int argc;
  495.     char *argv[];
  496.     {
  497.         strcpy(filename,"input.dat");
  498.  
  499.         read_data();
  500.         factor_A( A[][][], f[], B[][][])#;
  501.         form_sums( C[][][], sumCW[][], sumCz[])#;
  502.         form_J( sumCW[][], P[][] )#;
  503.         form_b( b[], sumCz[], fqp[] )#;
  504.         dist_lu( P[][], RowPerm[])#;
  505.         dist_solve( P[][], xqp[], b[], RowPerm[])#;
  506.         compute_xi(xqp[], x[])#;
  507.         print_results();
  508.     }
  509. }
  510.